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Abstract 

Evolutionary dynamics has been classically studied for homogeneous 
populations, but now there is a growing interest in the non-homogenous 
case. One of the most important models has been proposed in ffH. adapting 
to a weighted directed graph the process described in m- The Markov 
chain associated with the graph can be modified by erasing all non-trivial 
loops in its state space, obtaining the so-called Embedded Markov chain 
(EMC). The fixation probability remains unchanged, but the expected time 
to absorption (fixation or extinction) is reduced. In this paper, we shall 
use this idea to compute asymptotically the average hxation probability 
for complete bipartite graphs Kn,m- To this end, we firstly review some 
recent results on evolutionary dynamics on graphs trying to clarify some 
points. We also revisit the ‘Star Theorem’ proved in m for the star 
graphs Ki^rn- Theoretically, EMC techniques allow fast computation of 
the fixation probability, but in practice this is not always true. Thus, in 
the last part of the paper, we compare this algorithm with the standard 
Monte Carlo method for some kind of complex networks. 

Keywords: Evolutionary dynamics, Markov chain, Monte Carlo methods, fixa¬ 
tion probability, expected fixation time, star and bipartite graphs. 
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1 Introduction and motivation 

Population genetics studies the genetic composition of biological populations, 
and the changes in this composition that result from the action of four different 
processes: natural selection, random drift, mutation and migration. The modern 
evolutionary synthesis combines Darwin’s thesis on natural selection and Mendel’s 
theory of inheritance. According to this synthesis, the central object of study 
in evolutionary dynamics is the frequency distribution of the alternative forms 
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(allele) that a hereditary unit (gene) can take in a population evolving under 
these forces. 

Many mathematical models have been proposed to understand evolutionary 
process. Introduced in [12], the Moran model describes the change of gene 
frequency by random drift on a population of finite fixed size. This model has 
many variants, but we assume for simplicity that involved organisms are haploids 
with only two possible alleles a and A for a given locus. Suppose there is a 
single individual with a copy of the allele A. At each unit of time, one individual 
is chosen at random for reproduction and its clonal offspring replaces another 
individual chosen at random to die. To model natural selection, individuals with 
the advantageous allele A are assumed to have relative fitness r > 1 as compared 
with those with allele a of fitness 1. 

Evolutionary dynamics has been classically studied for homogeneous popula¬ 
tions, but it is a natural question to ask how non-homogeneous structures affect 
this dynamics. In im, a generalisation of the Moran process was introduced 
by arranging the population on a directed graph, see also |T3|, |IH| and [TU] . 
In this model, each vertex represents an individual in the population, and the 
offspring of each individual only replace direct successors, i.e. end-points of 
edges with origin in this vertex. The fitness of an individual represents again its 
reproductive rate which determines how often offspring takes over its neighbour 
vertices, although these vertices do not have to be replaced in an equiprobable 
way. The evolutionary process is described by the choice of stochastic matrix 
W = (wij) where Wij denotes the probability that individual i places its offspring 
into vertex j. In fact, further generalisations can be considered assuming that 
the probability above is proportional to the product of a weight Wij and the 
fitness of the individual i. In this case, W does not need to be stochastic, but 
non-negative. The fixation probability of the single individual i is the probability 
that the progeny of i takes over the whole population. Several interesting and 
important results are shown in m- 

• Different graph structures support different dynamical behaviours amplifying 
or suppressing the reproductive advantage of mutant individuals (with the 
advantageous allele A) over the resident individuals (with the disadvantageous 
allele a). 

• An evolutionary process on a weighted directed graph (G, W) is equivalent 

to a Moran process (i.e. there is a fixation probability well-defined for any 
individual, which coincides with the fixation probability in a homogeneous 
population) if and only if (G, W) is weight-balanced, i.e. for any vertex i the 
sum of the weights of entering edges W-{i) = J2j=i that of leaving 

edges w+{i) = equal. This is called the Circulation Theorem 

in m and m- 

As in the classical setting, mutant individuals will either become extinct or take 
over the whole population, reaching one of the two absorption states {extinction 
or fixation), when a finite population is arranged on an undirected graph or 
on a strongly connected directed graph (where two different vertices are always 
connected by an edge-path). Even in the first case, the fixation probability 
depends usually on the starting position of the mutant. The effect of this initial 
placement on mutant spread has been discussed in dij. 
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In the present paper, we start by summarising some fundamental ideas 
and results on evolutionary dynamics on graphs. In this context, most work 
involves computing the (average) fixation probability, but doing so in general 
requires solving a system of 2^ linear equations. In the example of the star 
graph described in m, like for other examples described in 0, m and [a, a 
high degree of symmetry reduces the size of the linear system to a set of 2N 
equations, which becomes asymptotically equivalent to a linear system with N 
equations. We revisit this example that will be useful in addressing the study of 
complete bipartite graphs. Another research direction has been to use Monte 
Carlo techniques to implement numerical simulations, but often limited to small 
graphs [1], small random modification of regular graphs m or graphs evolving 
under random drift m- 

Our aim is to show how to modify the stochastic process associated with a 
weighted directed graph to simplify the evolutionary process both analytically 
and numerically. Recall that an evolutionary process on a weighted directed 
graph (G, W) with N vertices is a Markov chain with 2^ states representing 
the vertex sets inhabited by mutant individuals and transition matrix P derived 
from W. The non-zero entries of P can be used to see the state space as a 
(weighted) directed graph. We call loop-erasing the loop suppression in this 
graph iS, avoiding to remain in the same state in two consecutive steps and 
providing the Embedded Markov chain (EMC) associated to the process. This 
technique is used here to compute asymptotically the average fixation probability 
for complete bipartite graphs, generalising the Star Theorem of [U, see also [J, 
[3] and [5T] . Expected time to absorption (fixation or extinction) of this EMC 
has been studied for circular, complete and star graphs in [7]. Here we compare 
numerically the expected absorption time of both chains on some kinds of 
complex networks. This method can be combined with other approximation 
methods (like the FPRAS method described in [5] for undirected graphs) to 
obtain a fast approximation scheme. 

The paper is organised as follows. In Section 2, we review the Moran model 
for homogeneous and non-homogeneous populations. In Section 3, we revisit the 
Star Theorem giving an alternative proof of it. In Section 4, we briefly explain 
the machinery of the loop-erasing method and we use this idea to describe the 
asymptotic behaviour of the fixation probability on the complete bipartite graphs 
family. At the end, in Section 5, we include some numerical experiments to 
evaluate the performance of the Monte Carlo method on both the standard and 
the loop-erased chains for different complex networks. 


2 Review of Moran process 

The Moran process models random drift and natural selection for finite homoge¬ 
neous populations m- As indicated before, we consider a haploid population 
of N individuals having only two possible alleles a and A for a given locus. At 
the beginning, all individuals have the allele a. Then one resident individual is 
chosen at random and replaced by a mutant having the neutral or advantageous 
allele A. At successive steps, one randomly chosen individual replicates with 
probability proportional to the fitness r > 1 and its offspring replaces one indi¬ 
vidual randomly chosen to be eliminated, see Figure Since the future state 
depends only on the present state, the Moran process is a Markov chain A„ 
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Figure 1: Classical Moran process 


with state space 5 = {0,..., N} representing the number of mutant individuals 
with the allele A at the time step n. This is a stationary process because the 
probability Pij = P[X„+i = j\Xn = i] to pass from i to j mutant individuals 
does not depend on the time n. In fact, the number of mutant individuals can 
change at most by one at each step and hence the transition matrix P = {Pij) 
is a tridiagonal matrix where Pij = 0 if j ^ i — 1, i + 1. As Po,o = Pn,n = 1, 
the states i = 0 and i = N are absorbing^ whereas the other states are transient. 

The fixation probability of i mutant individuals 

<^>^ = ^.(r) = F[3n > 0 : = N\Xo = i] 

is the solution of the system of linear equations: 

$0 = 0 

( 1 ) 

$^ = 1 


where = \ — Pi^i-i—Pi^i+i. In particular, the probability of a single mutant to 
reach fixation $i = $i(r) is usually referred to as the fixation probability in short. 
To solve Q, we define yt = $i—$i_i which verifies X]i=i Vi — $»—$o = 1- Then, 
dividing each side of Q by Pi,i+i, we have yi+i = where = Pi^i^i/Pi^i+i 
is the death-birth rate. It follows pi = $iri}=i7i) hence the fixation 
probability is 


$1 




( 2 ) 


See [ID], [ 22 ] and [H]. 

If neither of alleles a and A is advantageous reproductively, the random 
drift phenomenon is modelled by the Moran process with fitness r = 1, and (§ 
becomes $i = 1/A^. On the contrary, if mutant individuals with the allele A 
have fitness r ^ I according to the hypothesis of natural selection, then 7 ^ = 1/r 
and therefore 


$1 


N-l 




1 — r ^ 1 

- >1 - 

1 — r~^ r 


(3) 


Moran process on graphs 

The Moran process for non-homogenous populations represented by graphs was 
introduced in m- Like for finite homogenous populations, the first natural 
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question is to determine the chance that the offspring of a mutant individual 
having an advantageous allele spreads through the graph reaching any vertex. 
But this chance depends obviously on the initial position of the individual 
(see mm and the global graph structure may significantly modify the balance 
between random drift and natural selection observed in homogeneous populations 
as proved in m- 

Let G = {V,E) be a directed graph, where V = A^} is the set of 

vertices and E is the set of edges. We assume that G is finite, connected 
and simple graph (without loops or multiple edges). Thus i? is a subset of 
{(j, j) € V X V I i ^ j}. An evolutionary process on G is again a Markov chain, 
but each state is now described by a set of vertices S G S = V{V) inhabited by 
mutant individuals having a neutral or advantageous allele A. This reproductive 
advantage is measured by the fitness r > 1. The transition probabilities of 
this Markov chain are defined from a non-negative matrix W = (wij) satisfying 
Wij = 0 {i,j) ^ E. More precisely, the transition probability between two 

states S,S' G S is given by 



if \ 5 = b'}. 


SiGV\S ^b 

if^\^' = {j}. 

^ SjGV ^b + J2iev\S '“^b 

^E..,GS«'b+E*,,Gns^b 

if S' = S', 

^ SjGV ^b + J2iev\S '“^b 

0 

otherwise, 


where rJ 2 iesJ 2 jev''^i 3 + J2iev\sJ2jev reproductive 

weights of the mutant and resident individuals, equal to rliSI -|- — [S'! = 

A^ -I- (r — l)|iS'| when the matrix W is stochastic. Note that S is the vertex set 
of a directed graph Q where two states S and S' are joined by an edge if and 
only if Ps^s' 7^ 0. Thus, the Moran process on a weighted directed graph (G, W) 
is the random walk on Q defined by the 2^ x 2^ stochastic matrix P = {Ps,s')- 
The fixation probability of any set S inhabited by mutant individuals 


= $s(G, W,r) = P[3n > 0 : = I/IXq = S] 


is still obtained as the solution of the linear equation 


= d>, (5) 

which is analogous to 0 for the classical Moran process. As in this case. S' = 0 
and S = V are absorbing states, but there may be other states of this type, 
as well as other recurrent states, so the probability that resident or mutant 
individuals reach fixation can be strictly less than 1. However, it is well-known 
(see O Sec. 111.7]) that ([^ has a unique solution if the only recurrent states 
are 0 and V. Thus, the population will still reach one of the two absorbing 
states: extinction or fixation of mutant individuals. If there are other recurrent 
states, absorbing or not, ([^ will have further solutions if no other restrictions 
are imposed, see the two-sources digraph below. But the probability of reaching 
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Figure 2: Complete graph 
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Figure 3: Biased random walk 


V from those states is 0, so adding these boundary conditions, the uniqueness of 
the fixation probability remain true. 

In this context, the fixation probability depends on the starting position of 
the mutant in the graph. This justifies the following definition: for any weighted 
directed graph (G, W), we call average fixation probability the average 

1 ^ 

= $A(G,IF,r) = 

i=l 

Complete graph Let be the complete graph with vertex set V = 
{!,... ,N} and edge set E = {{i,j) G V x V \ i j}. The classical Moran 
process is the Moran process on G = defined by the stochastic matrix 
W = [wij) where Wij = if i j, see Figure 2 Since G is symmetric 

(i.e. the automorphism group Aut(G) acts transitive^ on V and E) and W is 
preserved by the action of Aut(G), <i>{q = for all i ^ j, and then 
for all i. 

Weight-balanced graph Assume that (G, W) is weight-balanced so that the 
sum of the weights of entering edges W-{i) = J2j=i that of leaving edges 

w_|_(f) = equal for any vertex i G V. According to the Circulation 

Theorem of m, the number of elements of each state of the Moran process 
on (G, W) ‘performs’ a biased random walk on the integer interval [0, A^] with 
forward bias r > 1 and absorbing states 0 and N^ see Figurej^ Reciprocally, if the 
Moran process on (G, W) reduces to this process, then (G, W) is weight-balanced. 


Two-sources digraph Let G be a directed graph consisting of two vertices 
(labelled 1 and 2) having leaving degree 1 and one vertex (labelled 3) having 
entering degree 2, see Figure [4(^ There are four recurrent states {!}, {2}, {1, 3} 
and {2, 3}, the average extinction probability is equal to 1/3, and the average 
fixation probability is equal to 0. Nonetheless, there is another state {1,2} 
having fixation probability equal to 1, see Figure [4(b)| 
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Figure 4: Two-sources digraph and its state space 



Figure 5: Star graph 


3 Star graphs revisited 

Lieberman et al. showed in m there are some graph structures, for example 
star structures, acting as evolutionary amplifiers favouring advantageous alleles. 
The evolutionary dynamics on stars graphs has been also studied in [3]. We 
revisit here this example that is useful to understand the role of symmetry for 
computing fixation probabilities. A star graph G consists of V = m -I-1 vertices 
labelled 0, 1 ,... ,m where only the centre 0 is connected with the peripheral 
vertices 1,..., m, see Figure Since Aut(G) acts transitively on the peripheral 
vertices, the state space reduces to a set of 2N ordered pairs. The fixation 
probability of the state (i, e) is denoted by 

=P[3n > 0 : = (m, l)|Vo = (z,e)], 

where i is the number of peripheral vertices inhabited by mutant individuals and 
£ € {0,1} indicates whether or not there is a mutant individual at the centre. 
The evolutionary dynamics of a star structure is described by the system of 
linear equations 

‘ho.o = 0 

d>..i = $,+1.1 + i^7i$,.o + (1 - P+, - 

$ 2.0 = P^o^i,l + Pip^i-1,0 + (1 - P+o ~ 

$m,l — 1 

since transitions exist only between state {i, 1) (resp. (i,0)) and states (i -I-1,1), 
(i,0) and (7 1) ioT i < m (resp. {i — 1,0), (i, 1) and (i,0) for i > 0), see Figurej^ 


( 6 ) 

(7) 


7 








0 0 


0 


0 0 0 


( 0 , 1 ) 


( 1 . 1 ) 


(i-1,1) 


(i,l) (i+1,1) /'\ (m-1,1) 


(m,l) 


( 0 . 0 ) 


(1.0) 


(i-1.0) (i.O) 


(j+1,0) \/ (m-l,0) 


(m,0) 


0 0 0 0 0 0 

Figure 6: State space of a star graph 


The non-trivial entries of P are given by 
P+, = P[X„+i = (i + 1,1)|X„ = (i, 1)] = 
O:i=P[X„+i = (i,0)|^n = (bl)] = 


m — I 


r{% -\- 1) -\- m — i m 
m — i 


r{i + 1) + m — i 


= P[X„+1 = (i, 1)|X„ = (*, 0)] = 


= P[X„+1 = (i -1,0) = (bO)] = - 


ri + m + 1 — i 

1 i 


ri + m + 1 — i m 


and 


1 _ P+ _ P- = + ^ 

m r(i + 1) + m — i 

1 P+ _ p- = + ^ 

m ri + m + 1 — i 


In particular, we have: 


‘&0.1 = 


-$ 


1.1 


r + TO 

Thus, the death/birth rates are given by 

P, 


and $i_o = 


rm + 1 


-$ 


1 . 1 - 


i.i w j -O.o 1 

7 i.i = = — and 7^,0 = -^ = —. 

T>+ r P+o rm 

Like for Q, the linear equations ([^ and Q reduce to 

TTl 

^i+1,1 ~ ‘hyi = 7 i,l(‘^’i,l ~ ‘^’i.o) = —(‘J’i.l ~ ‘J’i.o); 

r 

— ‘hyo = 77o(‘f’i,o — ‘I’i-i.o) = —(‘hyo — ‘hi-i.o)- 

rm 


From (10), it is easy to obtain the following identity: 

vO /In / rm \ »-i+i 

= > — {——t) ^ 


1=1 


j.i> 


( 8 ) 


(9) 

( 10 ) 


Vz = 1, . . . , TO. 


( 11 ) 


















Now, using ([^ and we have the following equation: 


$ 


i+1,1 


- = — 


«'*.l - 


rm 

rm + 1 rm 


1 / rm \ • 

m \rm + 1 / 




i-1,1 


l — Z 

n 

i=i 

m 


J_y-j 

rmJ 


, i-j + l 


r{rm + 1) 

i-2 

-E 

i=i 


$i.i - 


m / 1 \ 
r \rmJ 


/ rm \ 

\ rm + 1 / 

\rm + 1 / 


id 


where 


i-2 


lim y 

m—>-+oo 


i=i 


(^ V”' ( 

f rm \ 

\rm/ 

\rm + 1 / 


j-i+i 


1 = 0 . 


Thus, when m —i +oo, the peripheral process with fixation probabilities 
becomes more and more close to the Moran process determined by 


1 


$*+1,1 - $*,i = ^($*.1 - $*-i.i)- 

According to (§, the average fixation probability is 

* 1* /I ^ 

$A —-H-— (-: 

Vto d 


( 12 ) 


m 


rm 


^1, 


m+1 ’ Vto + 1 r + m m + 1 rm + 1/ 

+ 00 , becomes more and more close to the fixation 


m + 1 

and therefore as m 

probability of the Moran process determined by (12) having fitness > 1. In 


short, the star structure is a quadratic amplifier of selection m in the sense 
that the average fixation probability of a mutant individual with fitness r > 1 
converges to 

1 — 

which is the fixation probability of a mutant with fitness > 1 in the Moran 
process. We will say these two evolutionary processes are asymptotically equiva¬ 
lent. 


4 Loop-erasing on complete bipartite graphs 

Let us consider a Moran process on a weighted directed graph (G, W). This is a 
random walk on the directed graph Q whose vertex set is S and whose transition 
matrix P = {Ps,S') is given by ([^. Two states S,S' G S are connected by an 
edge in Q if and only if Ps.s' 7^ 0- Let G be the directed graph obtained by 
suppressing any loop in Q that connects a non-absorbing state S to itself. For 
any pair S, S' G S such that S is non-absorbing, the transition probability Ps,S' 
is replaced by 

^iS'\S = {j}orS\S' = {i}, 

Ps,S' = <^ 1 - ( 13 ) 

I 0 otherwise. 
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Figure 7: Loop-erasing method 


where 


TTs = Ps,s = 1 - ( X! ^S,S\{i} 

jev\s i&s 


(14) 


is the probability of staying one time in the state S. Equivalently to (13), 

Ps,S' = ^sPs,S' 


(15) 


n>0 


where the n-th power of its is the probability of staying n times in the 
state S. We say the random walk on the directed graph Q defined by the 
transition matrix P is obtained by loop-erasing from the Moran process on 
{G,W), see Figure]^ This is the Embedded Markov chain (EMC) with state 
space S obtained by forcing the Moran process on (G, W) to change of state in 
each step. The fixation probability of any set S inhabited by mutant individuals 
remains unchanged = $s, because the system of linear equations 

‘J’S = Ps,S^S + X! ^S,SU{j}^Su{j} + X! ^S,S\{i}^S\{t} 

jev\s ies 


can be rewritten as 


Ps,su{j}'^su{j}+ '^Ps,s\{i}^s\{i}■ 

Jev\s ies 

The biased random walk described in Figure arises by loop-erasing in any 
process equivalent to the Moran process. 

Assuming that 0 and V are the only recurrent states in Q, we know the 
population will reach one of these two absorbing states, fixation or extinction, 
from any other subset S G V inhabited by mutant individuals. Moreover, the 
transition matrix P admits a box decomposition 


/1 

0 

0\ 

b 

Q 

C 

Vo 

0 

1/ 


(16) 


For this type of absorbing Markov chain, the expected absorption time (i.e. the 
expected number of steps needed to go from the state S to one of the absorbing 
states 0 or F) is given by the system of linear equations 


7-5=^ Ps,S'TS' + 1 

S'gSt 


(17) 


10 










9 



(a) Complete bipartite graph ii' 2,3 

Figure 8: Complete bipartite graphs 



where St is the set of transient states, that is, different from 
the box decomposition (16), the equation 0 reduces to 


r = (/d-Q)-il = ^Qn, 


and V. Using 


(18) 


n>0 


where {Id — Q)~^ is the fundamental matrix of the Markov chain and 1 is the 
vector with all the coordinates equal to 1. We have similar identities for the 
Markov chain obtained by loop-erasing. Thus, using the obvious notation, the 
new expected absorption time is given by 

f = (19) 

n>0 


where {Id - Q) ^ = J2n>o Q". The vector ts represents the expected number 
of state transitions until absorption when the Moran process starts from a set 
S. This quantity has been studied in [7] for circular, complete and star graphs. 
Since transition may not happen at every step, the following result is clear: 

Proposition 4.1. Let t be the expected absorption time for the Moran process 
on a weighted directed graph {G, W). Let t be the expected absorption time for 
the process obtained by applying the loop-erasing method. Then for each transient 
state S G St, we have ts < ts- 

For unweighted and undirected graphs, Diaz et al. show in that, with high 
probability, the expected absorption time is bounded by a polynomial in N of 
order 3, 4 and 6 when r < 1, r > 1 and r = 1. They have also constructed a fully 
polynomial randomised approximation scheme for the probability of fixation 
and extinction. The loop-erasing method can be used to reduce the expected 
absorption time making the approximation of the fixation probability faster. We 
explore this path in Section 

Complete bipartite graph 

Now, we use the loop-erasing method to calculate the asymptotic fixation 
probability of any complete bipartite graph. Recall that a complete bipartite 
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graph is a graph whose vertices can be divided into two disjoint sets 

Vi = TTii} and V 2 = m 2 } such that every edge connects a vertex 

in Vi to one in V 2 . In particular, a star graph is a bipartite graph Km,i- The 
fixation probability for these graphs has been also studied in [9] and [21] . 

According to the Circulation Theorem, as any vertex has the same number 
of connections, the evolutionary process on the complete bipartite graph ATm.m 
is equivalent to the Moran process, so they have the same fixation probability, 
see Figures 8 (b) and |10(b)| 

For a bipartite graph Kmi,m 2 with mi 7 ^ m 2 , like for star graphs, the state 
space reduces to the product S = {0,1,..., mi} x {0,1,..., m 2 } where each 
ordered pair (i, j) G S indicates that there are i vertices in Vi and j vertices in 
V 2 inhabited by mutant individuals. The evolutionary dynamics is described by 
the system of linear equations 




under the boundary conditions $ 0,0 = 0 and = Ij where the transition 

probabilities are given by 


= P[^n+l = (* + l,j)\Xn = {i,j)] 
=P[X„+i = (*-l,j)|A„ = (z,i)] 
/^F=P[X„+i = (z,j + l)|A„ = (z,i)] 
=P[X„+i = (z,j-l)|A„ = (z,j)] 


rj 

mi — % 

r{i + j) + N -{i+j) 

mi 

m 2 - j 

i 

r{i+j) + N - (i+j) 

mi 

ri 

m 2 - j 

r{i+j) + N - {i+j) 

m 2 

mi — i 

j 

r{i + j) + N - {i + j) 

m 2 


and N = mi +m 2 - The subscript (z,j) denote the initial state, while the arrows 
and are guidelines indicating the the direction of corresponding edge 
for the directed graph structure on the state space (so that the next state is 
(* + l,j), (*-l,j), (bJ + l), or (z,j-l) respectively), see Figurej^ By applying 
the loop-erasing method, we obtain the following new transition probabilities: 


= 

i,m2 


rm2 


mi + rm2 


and 


pi _ 
^ 1,7712 


mi 


mi -I- rm2 


for the state (z,m 2 ) and 




P.F 


pi 


r(mi — i)jm2 

(mi -I- rm2){mi — i)j + z(m2 — j){rmi + m2) 

_ _ zm 2 (m 2 - j) __ 

(mi -I- rm2)(mi — i)j + z(m2 — j){rmi + m2) 
rzmi(m 2 — j) 

(mi -I- rm2)(mi — i)j + z(m2 — j){rmi + m2) 
_mi (mi - i)j _ 

(mi -I- rm2)(mi — i)j + z(m2 — j){rmi + m2) 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 


for the state (z,j) with 0 < j < m 2 . Like for star graphs, all the symmetries 
in complete bipartite graphs has been used to reduces the state space of the 
evolutionary process to the vertex set of the directed graph Q described in 
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Figure 9: State space of a bipartite graph 


Figure]^ But neither this reduced process (random walk on ^), nor the process 
obtained by loop-erasing (random walk on the directed graph Q obtained by 
suppressing every loop connecting a non-absorbing state with itself) admit global 
symmetries. Note that states ( 1 , 7712 ) in the last row are never related to states 
(i,j) with 1 < j < m 2 by automorphisms of the weighted directed graphs Q 
and Q. Nevertheless, we prove that the map sending the state (i, j) to the state 
{i,j — 1) with 1 < j < m 2 becomes more and more close to a symmetry of the 
Embedded Markov chain when mi —> -Poo. Using the previous calculation of 
the asymptotic fixation probability for a star graph, we obtain the following 
theorem, that is also illustrated numerically in Figure [T0| 

Theorem 4.2. Let $. 4 (Elmi,m 2 ; ^6 ^^6 average fixation probability of a single 

mutant individual having a neutral or advantageous allele A with fitness r > 1 
in a Moran process on a eomplete bipartite graph Kmi ,m 2 • Then 

lim 4>2i(-?^mi,m2U) = lim = 1 - ^ 

mi —>-+00 m—>-+oo 

where is the average fixation probability of a single mutant individual 

having a neutral or advantageous allele A with fitness r^ > 1 in the classical 
Moran process on Km ■ 


Proof. We start by observing that, according to (201, we have: 



rm2 

mi + rm2 
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(a) Average fixation probabilities for (b) Average fixation probabilities for 

^ 2 . 10 i ^2,50i ^ 2,100 and A' 2,200 ^lO.lOi ^10,50: ^10,100 and A' 10,200 


Figure 10: Average fixation probabilities for Moran processes on some bipartite 
graphs obtained from Monte Carlo methods. The same figures can be obtained 
from the loop-erasing method. 


for 1 < j < m2 and 


^ _ r(mi - i)jm2 _ 

_|_ rm 2 ){mi — i)j + i{m 2 — j)irmi + m2) 
r{mi — i){j — l)m 2 


(mi -I- rm2)(mi — i){j — 1 ) -f i{m2 — j + l)(rmi -|- m2) 
rm2 

i m2 — j , . 

mi -I- rm2 H-•-(rmi -I- m2) 

mi — i J 

rm2 

i m2— j + l , . 

mi -I- rm2 H-: •-^;-(rmi -I- m2) 


mi — i J — 1 

z / m2 — j -I- 1 m2— j 


< rm2 


= rm2 • 


< rm2 


-A 


mi — i \ J — 1 


J 


{rmi + m 2 ) 


m 2 


z»i - i j{j - 1 ) 


(mi -I- rm 2 )^ 
(rmi + m2) 


mi — i 


(mi -I- rm 2 )^ 

- m2 {rmi + m2) 


{mi + rm2y 


for 2 < j < m2 and for 1 < z < mi — 1. Assuming mi > 2 i, we have < 1 
and hence 


We deduce 


mi — i 


- m2(rmi + m 2 ) 


(mi + rm2)^ 


< rm2 


m2{rmi + m 2 ) 
{mi + rm2)^ 


lim = 0 


mi—>-+oo 


bi-i 


(24) 


for 2 < j < m2 and for z > 1. Similarly, using (211, we have = 0 for 
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1 < j < m 2 , P^rno = 0 for 1 < z < mi — 1, and 


_ _ zm2(m2 - j) _ 

(mi + rm2)(mi — i)j + i(m2 — j)irmi + m 2 ) 
_ m 2 

, . mi — i j 

(mi + rm2) -^- • -: + rmi + m 2 


m 2 -J 


< 


m 2 


rmi + m 2 

for 1 < j < mi — 1. As before, it follows: 


= Km2 = 0 

mi —>-+00 ’ 


(25) 


for 1 < j < m 2 and for each i > 1. Next, using (22 1 , we have PJ^- = 0 for 

rimi{m2 — j) 


1 < j < m 2 , Pi m 2 = 0 for 1 < z < mi — 1 , and 


P^ - = - 

(mi + rm 2 )(mi — i)j + z(m 2 — j){rmi + m 2 ) 

ri 


(mi + rm 2 ) 


mi — i 



^ + —{rmi + m 2 ) 


(mi + rm2) 


mi m 2 — J 


for 1 < 7 < mi — 1. Since i < ^ and — < —^^ when mi > 2z and 

— •J 2 rT?i m.o 777,0 — 7 — 

j < m 2 — 1 , we have 


PP < 


mi m 2 m2—J 

2m2ri 


and therefore 


(mi + rm 2 ) 


lim P/ = 0 

mi—>-+oo 


(26) 


for 1 < J < m 2 and for z > 1. Finally, we have: 


hm _ PP - PP_, = hm _ - P-.i +PCi- P^-i + Ph - Pl-i = 0 


bi-i 


from (24), (25) and (26). Arguing inductively on the integer z > 1, this implies 


that the Moran process on the bipartite graph Kmi,m 2 reduces asymptotically 
to the Moran process on the star ATmiP when mi —>■ +cx), and hence 

lim ^A{Krm,Tn 2 ,r) - $a(?’, Almi,l, ?") = 0, 

mi—>-+oo 


that proves the theorem. 


□ 


5 Numerical experiments in complex networks 

Proposition |4.1| says that the expected number of steps until absorption in the 
loop-erased Markov chain is smaller or equal than that in the standard one. 
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Figure 11: Average computation times (in seconds) for Moran processes on 
small-world networks (Watts-Strogatz /3-model) using EMC and SMC methods 
for all r going from 0 to 10 with step size of 0.25. 


At first glance, this seems to imply that Monte Carlo method on the EMC 
(EMC method from now on) will stop before Monte Carlo on the standard chain, 
(Standard Monte Carlo or SMC method from now on), but there is a subtle 
difference between what the method does theoretically and what the computer 
actually does. 

First of all, we need to construct a weighted directed graph Q having 2” states. 
It is almost always unfeasible when n is relatively large, but it is easier for highly 
symmetric graphs. So simulations reproduce how individuals randomly spawn 
and die. In the SMC method, at each step, the chance of selecting a mutant 
individual for reproduction is proportional to the fitness r > 1. This uniformity 
allows us to update the new transition probabilities in constant time. However, 
in the EMC method the probability of choosing each individual for reproduction 
depends not only on its fitness but also on the fitness of its neighbours. More 
precisely, the probability that a particular mutant individual v leaves offspring 
at a particular time is proportional to the number of resident neighbours of v 
at that time. Similarly, if r; is a resident individual, the probability of choosing 
it for reproduction is proportional to the number of mutant neighbours. Thus, 
if w is the neighbour of v chosen to die, the EMC method needs to update the 
transition probabilities of each neighbour of w. On some graphs, this may lead 
to longer computation times. 

We compared the amount of time it takes to end the simulations for the 
two methods in a series of well-known complex network models. All simulations 
were done on a computer running MacOS X 10.9.3 with a quad-core i5 at 
2.5GHz and 8Gb of RAM. Graph construction and manipulation was done in 
Sage/NetworkX jSJUO]; but the simulation routines were written in C. 

Small-world networks 

Small-world networks were introduced in |24j as a family of random graphs with 
some properties of real networks. The construction is as follows: consider a 
circular graph of order n and connect the k nearest neighbours. Now, each edge 
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Figure 12: Average computation times (in seconds) for Moran processes on 
scale-free networks (Barabasi-Albert preferential attachment model) using EMC 
and SMC methods for r going from 0 to 10 with step size of 0.25. 


uv in the previous graph can be replaced with another edge uw with probability 
p. The resulting graph may be disconnected. 

We did the following experiment to test the speed of the two methods: Fixed 
k = A, for any n € {10,20,..., 100} and p S {0, 0.05,...,!}, 

• we construct 10 random graphs with parameters n and p; 

• for each of these graphs, we compute the average fixation probability using 
both methods 3 times with 1000 trials for every fitness r varying from 0 to 10 
with step size of 0.25. 

Averaging the 30 running times of each method we get an average computation 
time for both algorithms on the family of small-world networks with the pre¬ 
scribed parameters. The results are shown in Figure 0 As can be seen, the 
EMC method performs better than the SMC method on this family of networks. 

Scale-free networks 

The previous family of graphs lacked a fairly common property of real networks, 
namely a power-law degree distribution. In [2] the preferential attachment model 
was developed to solve the shortcomings of previous models. Start with m 
vertices connected in no way. Then one single vertex is added and connected to 
the initial vertices to obtain a star. At successive steps, another single vertex is 
added and connected to m of the previous vertices with probability ‘proportional’ 
to the degree. After n — m steps, the graph has n vertices and (n — m)m edges. 
In the real world, one expect to have small m compared to n as the global 
population is large and people known just a very small portion of the population. 


We ran a similar experiment as for small-world networks, using both methods 
3 times for 10 random graphs with 1000 trials for every fitness r varying from 0 to 
10 in 40 evenly disposed steps, but new relevant parameters are now the order of 
the graph n G {100, 200,..., 1000} and the ratio m/n G {0.005, 0.01,..., 0.05}. 
Thus, a point {n,m/n) in the plot corresponds to the average computation time 
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Graph 

EMC time 

SMC time 

SMC/EMC 

Aipo 

0.98 

0.97 

0.99 

Ai,50 

23.79 

30.04 

1.26 

All,100 

181.48 

190.43 

1.05 

All,200 

1471.91 

1369.67 

0.93 

Ar2,io 

0.81 

0.88 

1.09 

K2,50 

12.52 

16.17 

1.29 

K2,100 

95.10 

98.49 

1.04 

K2,200 

753.86 

694.88 

0.92 

Afio.io 

0.91 

0.90 

0.99 

Afio.50 

4.72 

4.53 

0.96 

ATio.ioo 

29.68 

23.68 

0.80 

Ario ,200 

217.40 

154.32 

0.71 


Table 1: Average computation times (in seconds) for Moran processes on the 
star graphs of orders 11, 51, 101 and 201, and on the complete bipartite graphs 
K 2 ,W: Ar 2 ,ioOj ^ 2,2001 Aiopoj Aig^so, Aigpoo and Aio, 20 o- All simulations 

with 1000 trials and r going from 0 to 10 with step size of 0.25. The last column 
shows how many times faster EMC is than the SMC method. 


on the random family with parameters n and m = [n ■ m/nj. The size of the 
population n has been multiplied by 10 with respect to the size of the small-world 
networks in the previous sample because we need to consider a population with 
n < 100 individuals if m/n = 0.01 and n > 200 individuals if m/n = 0.005. The 
result of the simulation can be seen in Figure]^ The SMC method performance 
is specially bad on graphs obtained by Barabasi-Albert preferential attachment 
with m = 1, whereas EMC method performs badly on the models with ‘large’ 
m where high degree vertices appear. Star graphs can be interpreted as graphs 
obtained by Barabasi-Albert preferential attachment in a single step. According 
to this interpretation, SMC method should improve the average computation 
times obtained by the EMC method. In Table [l] we compare these times on 
some star and complete bipartite graphs with the same number of trials and 
values of the fitness r. 


Hierarchical networks 


In [15j , a deterministic network was introduced as a heuristic model of metabolic 


networks (it can be seen in Figure 13(a) I. The graph has a power law distribution 


of the degree (scale-free topology) and mean clustering coefficient non-decreasing 
with size. 

The network is constructed inductively. In the first step, we define the 
network i?o as the complete graph of order four K^. Fix one vertex as the central 
vertex. The rest of the vertices are external vertices. Now, take three copies of 
i?o and join their external vertices with the central one of the original i?o and 
their central vertices together making a big triangle. The central vertex of the 
resulting graph Ri is the central vertex of the original Rq. The external vertices 
of Ri are the vertices of the other copies of Rq. You can repeat the process as 
many times as needed, obtaining graphs Rn of order 
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- Moran process n—>^ + oo 

1—1 Ravasz Rq 

1—1 Ravasz i?j 

1—1 Ravasz i?2 

Ravasz 

— Ravasz R.^ 


0 2 4 6 8 10 



Figure 13: (a) Third step in the construction of a hierarchical network by [15], 
and (b) fixation probability for different construction steps compared with the 
Moran process as n —)■ +oo. Here 100 000 trials was carried out to obtain a less 
noisy approximation of the fixation probability. 


Table 2: Average computation times (in seconds) for Moran processes on the 
hierarchical model by m using EMC and SMC methods, both with 1000 trials 
and r going from 0 to 10 with step size of 0.25. The last column shows how 
many times faster EMC is than the SMC method. 


Step 

Graph order 

EMC time 

SMC time 

SMC/EMC 

0 

4 

0.76 

0.82 

1.07 

1 

16 

0.88 

0.94 

1.06 

2 

64 

1.92 

5.15 

2.68 

3 

256 

24.94 

141.24 

5.66 

4 

1024 

707.61 

6923.04 

9.78 


We computed the average fixation probability and the average computation 
time on this family for the same values of fitness and trials as before. The results 
can be seen in Eigure |13(b)| and Table As one can see, the EMC method 
outperforms the SMC method by an increasing factor on the order n of the 
network. 


6 Conclusion 

In this paper, we review some fundamental ideas and results on evolutionary 
dynamics as introduced in m generalising the classic process described in m- 
But we also give insights on one of the major problems in this theory, to estimate 
the average fixation probability of a mutant with relative fitness r on a given 
graph. Exact solutions have only been computed for a few families of graphs |3| 
as, generally, one should solve a linear equation systems of 2^ equations, where 
N is the order of the graph. Even asymptotic behaviour is tricky to compute. 

The erasing of loops in the state space is the geometrical counterpart of a well 
known device in Markov chains, which is the basis behind embedded processes 
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and which consists of forcing the processes to evolve in each iteration. It is 
rather obvious and well known that the expected fixation time is reduced and the 
fixation probability is unchanged by this procedure. In this paper, we use this 
idea to compute asymptotically the fixation probability for the class of complete 
bipartite graphs, generalising the result of m for the star graph. In this case, 
the high degree symmetry reduces the problem to a set of 2N equations, which 
is asymptotically equivalent to a simpler linear system of N equations. For 
complete bipartite graphs, after erasing all non-trivial loops, partial symmetries 
arise asymptotically in the Moran process and reduce the Moran process to the 
particular case of a star graph. This is an important step since it shed some new 
light on the asymptotic behaviour of the fixation on bipartite graphs, which has 
recently been dealt with from other points of view in and . 

In practice, the Monte Carlo on the Embedded Markov chain (EMC method) 
may need to make more computations than the Monte Carlo method on the 
standard chain (SMC method), as it needs to keep track of different probabilities 
(one per vertex) that should be computed at runtime. We tested the speed of the 
new method in some celebrated families of graphs: the small world networks |24j , 
preferential attachment networks [2] and hierarchical networks m- These tests 
show the EMC method defeats the SMC method on large families of graphs, but 
not in all examples, as transition probabilities on the loop-erased chain depends 
heavily on the actual state. At first, the appearance of high degree vertices 
might look like culprit for this problem, but this is not the case: the hierarchical 
network of m has extremely large degree on some vertices. Although it is still 
unknown what makes EMC method become slower, we believe this method could 
be applied successfully to real networks. 

Globally, we think the present paper represents substantial progress towards 
understanding the complexity behind evolutionary dynamics on graphs. 
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